Intro Datenanalyse 1

Jan-Philipp Kolb

27 April 2017

Warum R nutzen

Gründe für die Nutzung von R

Gründe

Übersicht - warum R

R lässt sich kombinieren…

R für SPSS Nutzer

install.packages("Rcmdr")
library("Rcmdr")

Bob Munich - R for SPSS and SAS Users

Die Popularität von R

R Nutzer rund um die Welt

R Welt

R Welt

Wo sind die aktivsten Nutzer?

Aktivität Nutzer

Aktivität Nutzer

Erwartungen und Anforderungen

Das kann diese Schulung vermitteln:

Erwartungen und Anforderungen II

Das kann sie nicht leisten:

R herunterladen:

http://www.r-project.org/

Probleme mit Excel

Weil andere Programme große Fehler haben:

Probleme mit Excel

Vergleich mit anderen Programmen

Dein Freund das GUI

Open Source Programm R

https://cran.r-project.org/

Graphisches User Interface

Aber die meisten Menschen nutzen einen Editor oder ein graphical user interface (GUI).

Aus den folgenden Gründen:

Verschiedene GUIs

rstudio

rstudio

Download der Unterlagen

Auf github sind alle Unterlagen für diesen Kurs zu finden.

Wie nutzt man github?

Rstudio

Aufgabe - Vorbereitung

date()
## [1] "Thu Apr 27 22:45:09 2017"
sessionInfo()
## R version 3.3.3 (2017-03-06)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 7 x64 (build 7601) Service Pack 1
## 
## locale:
## [1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252   
## [3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                   
## [5] LC_TIME=German_Germany.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] backports_1.0.5 magrittr_1.5    rprojroot_1.2   tools_3.3.3    
##  [5] htmltools_0.3.5 yaml_2.1.13     Rcpp_0.12.6     stringi_1.1.1  
##  [9] rmarkdown_1.4   knitr_1.15.1    stringr_1.2.0   digest_0.6.12  
## [13] evaluate_0.10

Grundlagen im Umgang mit der Sprache R

R ist eine Objekt-orientierte Sprache

Vektoren und Zuweisungen

b <- c(1,2) # erzeugt ein Objekt mit den Zahlen 1 und 2
mean(b) # berechnet den Mittelwert
## [1] 1.5

Mit den folgenden Funktionen können wir etwas über die Eigenschaften des Objekts lernen:

length(b) # b hat die Länge 2
## [1] 2

Objektstruktur

str(b) # b ist ein numerischer Vektor
##  num [1:2] 1 2

Funktionen im base-Paket

Funktion Bedeutung Beispiel
length() Länge length(b)
max() Maximum max(b)
min() Minimum min(b)
sd() Standardabweichung sd(b)
var() Varianz var(b)
mean() Mittelwert mean(b)
median() Median median(b)

Diese Funktionen brauchen nur ein Argument.

Funktionen mit mehr Argumenten

Andere Funktionen brauchen mehr:

Argument Bedeutung Beispiel
quantile() 90 % Quantile quantile(b,.9)
sample() Stichprobe ziehen sample(b,1)

Beispiel - Funktionen mit einem Argument

max(b)
## [1] 2
min(b)
## [1] 1
sd(b)
## [1] 0.7071068
var(b)
## [1] 0.5

Funktionen mit einem Argument

mean(b)
## [1] 1.5
median(b)
## [1] 1.5

Funktionen mit mehr Argumenten

quantile(b,.9)
## 90% 
## 1.9
sample(b,1) 
## [1] 1

Übersicht Befehle

http://cran.r-project.org/doc/manuals/R-intro.html

Aufgabe - Zuweisungen und Funktionen

Erzeugen Sie einen Vektor b mit den Zahlen von 1 bis 5 und berechnen Sie…

  1. den Mittelwert

  2. die Varianz

  3. die Standardabweichung

  4. die quadratische Wurzel aus dem Mittelwert

Datentypen und Indizieren

Verschiedene Datentypen

Datentyp Beschreibung Beispiel
numeric ganze und reele Zahlen 5, 3.462
logical logische Werte FALSE, TRUE
character Buchstaben und Zeichenfolgen “Hallo”

Quelle: R. Münnich und M. Knobelspieß (2007): Einführung in das statistische Programmpaket R

Verschiedene Datentypen

b <- c(1,2) # numeric
log <- c(T,F) # logical
char <-c("A","b") # character
fac <- as.factor(c(1,2)) # factor

Mit str() bekommt man den Objekttyp.

Indizieren eines Vektors:

A1 <- c(1,2,3,4)
A1
## [1] 1 2 3 4
A1[1]
## [1] 1
A1[4]
## [1] 4
A1[1:3]
## [1] 1 2 3
A1[-4]
## [1] 1 2 3

data.frames

Beispieldaten generieren:

AGE <- c(20,35,48,12)
SEX <- c("m","w","w","m")

Diese beiden Vektoren zu einem data.frame verbinden:

Daten <- data.frame(Alter=AGE,Geschlecht=SEX)

Anzahl der Zeilen/Spalten herausfinden

nrow(Daten) # Zeilen
## [1] 4
ncol(Daten) # Spalten
## [1] 2

Indizieren

Indizieren eines dataframe:

AA <- 4:1
A2 <- cbind(A1,AA)
A2[1,1]
## A1 
##  1
A2[2,]
## A1 AA 
##  2  3
A2[,1]
## [1] 1 2 3 4
A2[,1:2]
##      A1 AA
## [1,]  1  4
## [2,]  2  3
## [3,]  3  2
## [4,]  4  1

Matrizen und Arrays

A <- matrix(seq(1,100), nrow = 4)
dim(A)
## [1]  4 25

Ein Array erzeugen

A3 <- array(1:8,c(2,2,2))
A3
## , , 1
## 
##      [,1] [,2]
## [1,]    1    3
## [2,]    2    4
## 
## , , 2
## 
##      [,1] [,2]
## [1,]    5    7
## [2,]    6    8

Indizieren eines Array

A3[,,2]
##      [,1] [,2]
## [1,]    5    7
## [2,]    6    8

Listen

Indizieren einer Liste

A4 <- list(A1,1)
A4
## [[1]]
## [1] 1 2 3 4
## 
## [[2]]
## [1] 1
A4[[2]]
## [1] 1

Logische Operatoren

# Ist 1 größer als 2?
1>2
## [1] FALSE
1<2
## [1] TRUE
1==2
## [1] FALSE

Operatoren um Subset für Datensatz zu bekommen

Diese Operatoren eignen sich gut um Datensätze einzuschränken

Daten
##   Alter Geschlecht
## 1    20          m
## 2    35          w
## 3    48          w
## 4    12          m
Daten[AGE>20,]
##   Alter Geschlecht
## 2    35          w
## 3    48          w

Datensätze einschränken

Daten[SEX=="w",]
##   Alter Geschlecht
## 2    35          w
## 3    48          w
# gleiches Ergebnis:
Daten[SEX!="m",]
##   Alter Geschlecht
## 2    35          w
## 3    48          w

Weitere wichtige Optionen

# Ergebnis in ein Objekt speichern
subDat <- Daten[AGE>20,]
# mehrere Bedingeungen können mit
# & verknüpft werden:
Daten[AGE>18 & SEX=="w",]
##   Alter Geschlecht
## 2    35          w
## 3    48          w

Sequenzen

# Sequenz von 1 bis 10
1:10
##  [1]  1  2  3  4  5  6  7  8  9 10
Daten[1:3,]
##   Alter Geschlecht
## 1    20          m
## 2    35          w
## 3    48          w

Weitere Sequenzen

seq(-2,8,by=1.5)
## [1] -2.0 -0.5  1.0  2.5  4.0  5.5  7.0
a <-seq(3,12,length=12)

b <- seq(to=5,length=12,by=0.2)

d <- 1:10
d <- seq(1,10,1)
d <- seq(length=10,from=1,by=1)

Wiederholungen

# wiederhole 1 10 mal
rep(1,10)
##  [1] 1 1 1 1 1 1 1 1 1 1
rep("A",10)
##  [1] "A" "A" "A" "A" "A" "A" "A" "A" "A" "A"

Die Funktion paste

?paste
paste(1:4)
## [1] "1" "2" "3" "4"
paste("A", 1:6, sep = "")
## [1] "A1" "A2" "A3" "A4" "A5" "A6"

Wie bekommt man Hilfe?

Wie bekommt man Hilfe?

help.start()
help(name)
?mean
example(lm)

Nutzung Suchmaschinen

R-project + Was ich schon immer wissen wollte

Stackoverflow

Ein Schummelzettel - Cheatsheet

https://www.rstudio.com/resources/cheatsheets/

Modularer Aufbau von R

Modularer Aufbau

Modularer Aufbau

install.packages("lme4")

library(lme4)

Installation von Paketen mit RStudio

Vorhandene Pakete und Installation

Übersicht viele nützliche Pakete:

Die wichtigsten Pakete zur Visualisierung mit R:

Pakete Regression

Big Data

Weitere interessante Pakete

Paket für den Import/Export - foreign

Pakete von Github installieren

install.packages("devtools")
library(devtools)

install_github("hadley/ggplot2")

Wie bekomme ich einen Überblick

Aufgabe - Zusatzpakete

Gehen Sie auf <cran.r-project.org> und suchen Sie in dem Bereich, wo die Pakete vorgestellt werden, nach Paketen,…

Datenimport

Datenimport

Dateiformate in R

Formate - base package

R unterstützt von Haus aus schon einige wichtige Formate:

Der Arbeitsspeicher

So findet man heraus, in welchem Verzeichnis man sich gerade befindet

getwd()

So kann man das Arbeitsverzeichnis ändern:

Man erzeugt ein Objekt in dem man den Pfad abspeichert:

main.path <- "C:/" # Beispiel für Windows
main.path <- "/users/Name/" # Beispiel für Mac
main.path <- "/home/user/" # Beispiel für Linux

Und ändert dann den Pfad mit setwd()

setwd(main.path)

Bei Windows ist es wichtig Slashs anstelle von Backslashs zu verwenden.

Alternative - Arbeitsspeicher

Import von Excel-Daten

library(foreign)
?read.csv
?read.csv2

CSV Dateien einlesen

Zunächst muss das Arbeitsverzeichnis gesetzt werden, in dem sich die Daten befinden:

Dat <- read.csv("schuldaten_export.csv")

Wenn es sich um Deutsche Daten handelt:

Dat <- read.csv2("schuldaten_export.csv")

SPSS Dateien einlesen

Dateien können auch direkt aus dem Internet geladen werden:

link<- "http://www.statistik.at/web_de/static/
mz_2013_sds_-_datensatz_080469.sav"

?read.spss
Dat <- read.spss(link,to.data.frame=T)

stata Dateien einlesen

MZ02 <- read.dta("MZ02.dta")

Das Paket rio

install.packages("rio")
library("rio")
x <- import("mtcars.csv")
y <- import("mtcars.rds")
z <- import("mtcars.dta")

Datenmanagement ähnlich wie in SPSS oder Stata

install.packages("Rz")
library(Rz)

Weitere Alternative Rcmdr

install.packages("Rcmdr")
library(Rcmdr)

Aufgabe - Datenimport

Datenexport

Datenexport

R’s Exportformate

Beispieldatensatz erzeugen

A <- c(1,2,3,4)
B <- c("A","B","C","D")

mydata <- data.frame(A,B)

Überblick Daten Import/Export

save(mydata, file="mydata.RData")

Daten in Excel Format abspeichern

write.csv(mydata,file="mydata.csv") 
library(xlsx)
write.xlsx(mydata,file="mydata.xlsx")

Daten in stata Format abspeichern

library(foreign)
write.dta(mydata,file="mydata.dta")

Auch zum Export eignet sich das rio Paket

library("rio")

export(mtcars, "mtcars.csv")
export(mtcars, "mtcars.rds")
export(mtcars, "mtcars.dta")

Exkurs: Datenquellen

Datenzugang

Datenquellen

library("FAOSTAT")
library(survey)
data(nhanes)

Das R-Paket datasets

library(datasets)

Beispiel Erdbeben Datensatz:

head(quakes)

Datensatz zum US Zensus

library(UScensus2010)

Weltbank Daten

WDI - World Development Indicators (World Bank) - Einführung in das Paket

library(WDI)
WDIsearch('gdp')[1:10,]

Nutzung von WDI Daten

dat <-  WDI(indicator='NY.GDP.PCAP.KD', country=c('MX','CA','US'), start=1960, end=2012)
head(dat)

Erste Grafik mit WDI Daten

OpenStreetMap

OpenStreetMap (OSM) ist ein kollaboratives Projekt um eine editierbare Weltkarte zu erzeugen.

Wikipedia - OpenStreetMap

Download von OpenStreetMap Daten

library(osmar)
api <- osmsource_api()
library(ggmap)
cityC <- geocode("Berlin",source="google")
bb <- center_bbox(cityC$lon,cityC$lat,1000, 1000)
uaBerlin <- get_osm(bb, source = api)

TwittR

library(twitteR)
library(streamR)

http://www.r-bloggers.com/mapping-the-world-with-tweets-including-a-gif-without-cats-and-a-shiny-app/

worldHires Daten

library(mapdata)
data(worldHiresMapEnv)
map('worldHires', col=1:10)

Historische Daten

library(HistData)
data(Arbuthnot)

GDELT Daten

library(GDELTtools)
test.filter <- list(ActionGeo_ADM1Code=c("NI", "US"), ActionGeo_CountryCode="US")
test.results <- GetGDELT(start.date="1979-01-01", end.date="1979-12-31",
                         filter=test.filter)

Andere Datenquellen

link1 <- "http://openflights.svn.sourceforge.net/viewvc/openflights/
openflights/data/airports.dat"
airport <- read.csv(link1, header = F)

link2 <- "http://openflights.svn.sourceforge.net/viewvc/openflights/
openflights/data/routes.dat"
route <- read.csv(link2, header = F)
link <- "http://www.fa-technik.adfc.de/code/opengeodb/DE9.tab"
info <- read.csv(link,sep="\t",header=F)

Weitere Quellen

Datenanalyse

Streuungsmaße

Im base Paket sind die wichtigsten Streuungsmaße enthalten:

ab <- rnorm(100); var(ab)
## [1] 1.072547
sd(ab); range(ab)
## [1] 1.035638
## [1] -2.520023  2.572660

Extremwerte

min(ab)
## [1] -2.520023
max(ab)
## [1] 2.57266

Fehlende Werte

ab[10] <- NA

var(ab)
## [1] NA

Bei fehlenden Werten muss ein weiteres Argument mitgegeben werden:

var(ab,na.rm=T)
## [1] 1.078947

Häufigkeiten und gruppierte Kennwerte

x <- sample(1:10,100,replace=T)

table(x)
## x
##  1  2  3  4  5  6  7  8  9 10 
## 11  7 10 11  8 14  7  9 16  7

Tabellieren - weiteres Beispiel

musician <- sample(c("yes","no"),100,replace=T)
?table
table(x)
## x
##  1  2  3  4  5  6  7  8  9 10 
## 11  7 10 11  8 14  7  9 16  7
table(x,musician)
##     musician
## x    no yes
##   1   6   5
##   2   4   3
##   3   7   3
##   4   6   5
##   5   3   5
##   6   6   8
##   7   4   3
##   8   5   4
##   9   6  10
##   10  5   2

Eine weitere Tabelle

data(esoph)
table(esoph$agegp)
## 
## 25-34 35-44 45-54 55-64 65-74   75+ 
##    15    15    16    16    15    11

Häufigkeitstabellen

Die Funktion prop.table()

table(esoph$agegp,esoph$alcgp)
##        
##         0-39g/day 40-79 80-119 120+
##   25-34         4     4      3    4
##   35-44         4     4      4    3
##   45-54         4     4      4    4
##   55-64         4     4      4    4
##   65-74         4     3      4    4
##   75+           3     4      2    2

Die Funktion prop.table

?prop.table
prop.table(table(esoph$agegp,
esoph$alcgp),1)
##        
##         0-39g/day     40-79    80-119      120+
##   25-34 0.2666667 0.2666667 0.2000000 0.2666667
##   35-44 0.2666667 0.2666667 0.2666667 0.2000000
##   45-54 0.2500000 0.2500000 0.2500000 0.2500000
##   55-64 0.2500000 0.2500000 0.2500000 0.2500000
##   65-74 0.2666667 0.2000000 0.2666667 0.2666667
##   75+   0.2727273 0.3636364 0.1818182 0.1818182

Die aggregate Funktion

aggregate(state.x77,by=list(state.region),mean)
##         Group.1 Population   Income Illiteracy Life Exp    Murder  HS Grad
## 1     Northeast   5495.111 4570.222   1.000000 71.26444  4.722222 53.96667
## 2         South   4208.125 4011.938   1.737500 69.70625 10.581250 44.34375
## 3 North Central   4803.000 4611.083   0.700000 71.76667  5.275000 54.51667
## 4          West   2915.308 4702.615   1.023077 71.23462  7.215385 62.00000
##      Frost      Area
## 1 132.7778  18141.00
## 2  64.6250  54605.12
## 3 138.8333  62652.00
## 4 102.1538 134463.00

x: ein oder mehrere Beobachtungsvektor(en) für den der Kennwert berechnet werden soll

by: eine oder mehrere bedingende Variable(n)

FUN: die Funktion welche den Kennwert berechnet (z.B. mean oder sd)

Beispieldatensatz - apply Funktion

ApplyDat <- cbind(1:4,runif(4),rnorm(4))
apply(ApplyDat,1,mean)
## [1] 0.2174516 0.3037641 1.0847166 1.6806857
apply(ApplyDat,2,mean)
## [1]  2.5000000  0.4733878 -0.5084242

Die Funktion apply

apply(ApplyDat,1,var)
## [1] 0.6143683 2.4650734 3.4903401 4.0587141
apply(ApplyDat,1,sd)
## [1] 0.7838165 1.5700552 1.8682452 2.0146251
apply(ApplyDat,1,range)
##            [,1]      [,2]       [,3]      [,4]
## [1,] -0.5676269 -1.098575 -0.7326386 0.3651438
## [2,]  1.0000000  2.000000  3.0000000 4.0000000
apply(ApplyDat,1,length)
## [1] 3 3 3 3

Argumente der Funktion apply

Die Funktion tapply

ApplyDat <- data.frame(Income=rnorm(5,1400,200),
                       Sex=sample(c(1,2),5,replace=T))
ApplyDat
##     Income Sex
## 1 1336.306   2
## 2 1231.769   2
## 3 1374.520   1
## 4 1430.820   1
## 5 1280.671   1

Beispiel Funktion tapply

tapply(ApplyDat$Income,ApplyDat$Sex,mean)
##        1        2 
## 1362.004 1284.037
tapply(ApplyDat$Income,
       ApplyDat$Sex,function(x)x)
## $`1`
## [1] 1374.520 1430.820 1280.671
## 
## $`2`
## [1] 1336.306 1231.769

Aufgabe - Apply Funktion anwenden

Zurück zur Gliederung.

Einfache Grafiken

Ein Plot sagt mehr als 1000 Worte

Plot ist nicht gleich Plot

CRAN Task Views

install.packages("ctv")
library("ctv")
install.views("Bayesian")

Task View zu Thema Graphiken

Datensatz

library(mlmRev)
data(Chem97)

Histogramm - Die Funktion hist()

Wir erstellen ein Histogramm der Variable gcsescore:

?hist
hist(Chem97$gcsescore)

Graphik speichern

Befehl um Graphik zu speichern

png("Histogramm.png")
hist(Chem97$gcsescore)
dev.off()

Histogramme

Argument Bedeutung Beispiel
main Überschrift main=“Hallo Welt”
xlab x-Achsenbeschriftung xlab=“x-Werte”
ylab y-Achsenbeschriftung ylab=“y-Werte”
col Farbe col=“blue”

Histogramm

hist(Chem97$gcsescore,col="blue",
     main="Hallo Welt",ylab="y-Werte", xlab="x-Werte")

Weitere Argumente:

?plot
# oder
?par

Barplot

tabScore <- table(Chem97$score)
barplot(tabScore)

Barplots und barcharts

barplot(tabScore)

Mehr Farben:

barplot(tabScore,col=rgb(0,0,1))

Grüne Farbe

barplot(tabScore,col=rgb(0,1,0))

Rote Farbe

barplot(tabScore,col=rgb(1,0,0))

Transparent

barplot(tabScore,col=rgb(1,0,0,.3))

Boxplot

?boxplot

Horizontaler Boxplot

boxplot(Chem97$gcsescore,
horizontal=TRUE)

Gruppierte Boxplots

Beispiel grupierter Boxplot

boxplot(Chem97$gcsescore~Chem97$gender)

Alternativen zu Boxplot

Violinplot

# Beispieldaten erzeugen
x <- rnorm(100)
y <- rnorm(100)

Die Bibliothek vioplot

library(vioplot)
plot(x, y, xlim=c(-5,5), ylim=c(-5,5))
vioplot(x, col="tomato", horizontal=TRUE, at=-4, 
        add=TRUE,lty=2, rectCol="gray")
vioplot(y, col="cyan", horizontal=FALSE, at=-4, 
        add=TRUE,lty=2)

vioplot - Das Ergebnis

Alternativen zum Boxplot

library(beanplot)
par(mfrow = c(1,2))
boxplot(count~spray,data=InsectSprays,col="blue")
beanplot(count~spray,data=InsectSprays,col="orange")

Grafiken für bedingte, bi- und multivariate Verteilungen

Scatterplots

Aufgabe - einfache Grafiken

Zurück zur Gliederung.

Zusammenhang

Edgar Anderson’s Iris Daten

data(iris)
head(iris)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa

petal length and width - Blütenblatt Länge und Breite

sepal length and width - Kelchblatt Länge und Breite

Zusammenhang zwischen stetigen Variablen

# Pearson Korrelationskoeffizient
cor(iris$Sepal.Length,iris$Petal.Length)
## [1] 0.8717538

Zusammenhang zwischen mehreren Variablen

pairs(iris[,1:4])

Zusammenhang zwischen mehreren Variablen

library("psych")
pairs.panels(iris[1:4],bg=c("red","yellow","blue")
[iris$Species],pch=21,main="")

Verschiedene Korrelationskoeffizienten

# Pearson Korrelationskoeffizient
cor(iris[,1:4]) 
##              Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length    1.0000000  -0.1175698    0.8717538   0.8179411
## Sepal.Width    -0.1175698   1.0000000   -0.4284401  -0.3661259
## Petal.Length    0.8717538  -0.4284401    1.0000000   0.9628654
## Petal.Width     0.8179411  -0.3661259    0.9628654   1.0000000
# Kendall's tau (Rangkorrelation)
cor(iris[,1:4], method = "kendall") 
##              Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length   1.00000000 -0.07699679    0.7185159   0.6553086
## Sepal.Width   -0.07699679  1.00000000   -0.1859944  -0.1571257
## Petal.Length   0.71851593 -0.18599442    1.0000000   0.8068907
## Petal.Width    0.65530856 -0.15712566    0.8068907   1.0000000
# Spearman's rho (Rangkorrelation)
cor(iris[,1:4], method = "spearman") 
##              Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length    1.0000000  -0.1667777    0.8818981   0.8342888
## Sepal.Width    -0.1667777   1.0000000   -0.3096351  -0.2890317
## Petal.Length    0.8818981  -0.3096351    1.0000000   0.9376668
## Petal.Width     0.8342888  -0.2890317    0.9376668   1.0000000

Zusammenhang zwischen kategorialen Variablen

Levelplot

library("lattice")
library("AER")
data(BankWages)
levelplot(table(BankWages$education,BankWages$job))

Visualisierung von Zusammenhängen zwischen kategorialen Variablen

mosaicplot(~ Sex + Age + Survived, 
           data = Titanic, color = TRUE)

Shading

Flächen werden entsprechend der Residuen eingefärbt:

mosaicplot(~ Sex + Age + Survived, 
           data = Titanic, shade = TRUE)

Literatur zu Zusammenhangsmaßen

Sachs - Angewandte Statistik mit R

Das lattice Paket

Das lattice-Paket

It is designed to meet most typical graphics needs with minimal tuning, but can also be easily extended to handle most nonstandard requirements.

http://stat.ethz.ch/R-manual/R-devel/library/lattice/html/Lattice.html

Histogramm mit Lattice

library("lattice");library("mlmRev")
data(Chem97)
histogram(~ gcsescore, data = Chem97)

Histogramm mit Lattice

  histogram(~ gcsescore | factor(score),data = Chem97)

Die Dichte mit Lattice zeichnen

densityplot(~ gcsescore | factor(score), Chem97, 
    groups=gender,plot.points=FALSE,auto.key=TRUE)

Einführung in das Paket lattice

Boxplot mit Lattice zeichnen

bwplot(factor(score) ~ gcsescore | gender, Chem97)

Boxplot mit Lattice zeichnen

bwplot(gcsescore ~ gender | factor(score), Chem97,
 layout = c(6, 1))

Univariate Plots

barchart(yield ~ variety | site, data = barley,
         groups = year, layout = c(1,6), stack = TRUE,
         auto.key = list(space = "right"),
         ylab = "Barley Yield (bushels/acre)",
         scales = list(x = list(rot = 45)))

Densityplot

densityplot( ~ height | voice.part, data = singer, layout = c(2, 4),  
            xlab = "Height (inches)", bw = 5)

Bivariate Plots

qq(gender ~ gcsescore | factor(score), Chem97,
   f.value = ppoints(100), type = c("p", "g"), aspect = 1)

xyplot

xyplot(Sepal.Length + Sepal.Width ~ Petal.Length + Petal.Width | Species,
       data = iris, scales = "free", layout = c(2, 2),
       auto.key = list(x = .6, y = .7, corner = c(0, 0)))

Multivariate Plots

splom(~iris[1:4], groups = Species, data = iris)

super.sym <- trellis.par.get("superpose.symbol")
splom(~iris[1:4], groups = Species, data = iris,
      panel = panel.superpose,
      key = list(title = "Three Varieties of Iris",
                 columns = 3, 
                 points = list(pch = super.sym$pch[1:3],
                 col = super.sym$col[1:3]),
                 text = list(c("Setosa", "Versicolor", "Virginica"))))

parallelplot

parallelplot(~iris[1:4] | Species, iris)

Lattice Befehle

Aufgabe - OECD Datensatz

data <- read.csv("oecd.csv",header = TRUE)

Zurück zur Gliederung.

Die lineare Regression

Die lineare Regression

Maindonald - DataAnalysis

Lineare Regression in R - Beispieldatensatz

John H. Maindonald and W. John Braun

DAAG - Data Analysis and Graphics Data and Functions

install.packages("DAAG")
library("DAAG")
data(roller)

help on roller data:

?roller

Das lineare Regressionsmodell in R

Schätzen eines Regressionsmodells:

roller.lm <- lm(depression ~ weight, data = roller)

So bekommt man die Schätzwerte:

summary(roller.lm)
## 
## Call:
## lm(formula = depression ~ weight, data = roller)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.180 -5.580 -1.346  5.920  8.020 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  -2.0871     4.7543  -0.439  0.67227   
## weight        2.6667     0.7002   3.808  0.00518 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.735 on 8 degrees of freedom
## Multiple R-squared:  0.6445, Adjusted R-squared:  0.6001 
## F-statistic:  14.5 on 1 and 8 DF,  p-value: 0.005175

Falls das Modell ohne Intercept geschätzt werden soll:

lm(depression ~ -1 + weight, data = roller)
## 
## Call:
## lm(formula = depression ~ -1 + weight, data = roller)
## 
## Coefficients:
## weight  
##  2.392

Summary des Modells

summary(roller.lm)
## 
## Call:
## lm(formula = depression ~ weight, data = roller)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.180 -5.580 -1.346  5.920  8.020 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  -2.0871     4.7543  -0.439  0.67227   
## weight        2.6667     0.7002   3.808  0.00518 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.735 on 8 degrees of freedom
## Multiple R-squared:  0.6445, Adjusted R-squared:  0.6001 
## F-statistic:  14.5 on 1 and 8 DF,  p-value: 0.005175

R arbeitet mit Objekten

predict(roller.lm) # Vorhersage
##         1         2         3         4         5         6         7 
##  2.979669  6.179765  6.713114 10.713233 12.046606 14.180002 14.980026 
##         8         9        10 
## 18.180121 24.046962 30.980502
resid(roller.lm) # Residuen
##          1          2          3          4          5          6 
## -0.9796695 -5.1797646 -1.7131138 -5.7132327  7.9533944  5.8199976 
##          7          8          9         10 
##  8.0199738 -8.1801213  5.9530377 -5.9805017

Residuenplot

plot(roller.lm,1)

Residuenplot

plot(roller.lm,2)

Linkliste - lineare Regression

Aufgabe - lineare Regression

Mietspiegel München

Beschrieben wird Wegstrecke, dreier Spielzeugautos die in unterschiedlichen Winkeln Rampe herunterfuhren.

  1. Lesen Sie den Datensatz toycars in einen dataframe data ein und wandeln Sie die Variable car des Datensatzes in einen Faktor (as.factor) um.

  2. Erstellen Sie drei Boxplots, die die zurückgelegte Strecke getrennt nach dem Faktor car darstellen.

  3. Schätzen Sie für jedes der 3 Autos separat die Parameter des folgenden linearen Mo dells mit Hilfe der Funktion lm()

\[ distance_i= \beta_0 + \beta_1 \cdot angle_i + \epsilon_i\]

  1. Überprüfen Sie deskriptiv die Anpassung der drei Modelle, indem Sie die Regressionger- ade in einen Plot von distance gegen angle einfügen. Deutet das \[ R^2 \] jeweils auf eine gute Modellanpassung hin?

Zurück zur Gliederung.

Die logistische Regression

Agresti - Categorical Data Analysis (2002)

Faraway Bücher zu Regression in R

Binäre AVs mit glm

Beispieldaten für die logistische Regression

install.packages("HSAUR")
library("HSAUR")
data("plasma", package = "HSAUR")

Logistische Regression mit R

cdplot(ESR ~ fibrinogen, data = plasma)

Logistische Regression mit R

plasma_glm_1 <- glm(ESR ~ fibrinogen, data = plasma, 
                    family = binomial())

Weitere Beispieldaten

install.packages("faraway")
library("faraway")
data(orings)
temp damage
53 5
57 1
58 1

Generalisierte Regression mit R - weitere Funktionen

probitmod <- glm(cbind(damage,6-damage) ~ temp, 
    family=binomial(link=probit), orings)
modp <- glm(Species ~ .,family=poisson,gala)
library("MASS")
house.plr<-polr(Sat~Infl,weights=Freq,data=housing)

Linkliste - logistische Regression

Aufgabe - Datenanalyse

Zurück zur Gliederung.

Faktoren in R

Was sind Faktoren in R

http://www.stat.berkeley.edu/~s133/factors.html

Beispiel Definition von Faktoren

data <- c(1,2,2,3,1,2,3,3,1,2,3,3,1)
fdata <- factor(data)
fdata
##  [1] 1 2 2 3 1 2 3 3 1 2 3 3 1
## Levels: 1 2 3
rdata <- factor(data,labels=c("I","II","III"))
rdata
##  [1] I   II  II  III I   II  III III I   II  III III I  
## Levels: I II III

Weitere Möglichkeit der Definition

levels(fdata) <- c('I','II','III')
fdata
##  [1] I   II  II  III I   II  III III I   II  III III I  
## Levels: I II III

Beispiel Monate

mons <- c("March","April","January",
          "November","January","September",
          "October","September","November",
          "August","January","November",
          "November","February","May",
          "August","July","December",
          "August","August","September",
          "November","February","April")
mons <- factor(mons)
table(mons)
## mons
##     April    August  December  February   January      July     March 
##         2         4         1         2         3         1         1 
##       May  November   October September 
##         1         5         1         3

Beispiel: ordered factor

mons <- factor(mons,levels=c("January","February",
                             "March","April","May","June",
                             "July","August","September",  
                             "October","November",
                             "December"),
               ordered=TRUE)

mons[1] < mons[2]
## [1] TRUE

Rücktransformation

fert <- c(10,20,20,50,10,20,10,50,20)
fert <- factor(fert,levels=c(10,20,50),ordered=TRUE)
fert
## [1] 10 20 20 50 10 20 10 50 20
## Levels: 10 < 20 < 50
mean(fert)
## [1] NA
mean(as.numeric(levels(fert)[fert]))
## [1] 23.33333

Tabellen mit Faktoren

lets <- sample(letters,size=100,replace=TRUE)
lets <- factor(lets)
table(lets[1:5])
## 
## a b c d e f g h i j k l m n o p q r s t u v w x y z 
## 0 1 1 0 0 0 0 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0

Grafiken mit ggplot

Das Paket ggplot2

Basiseinführung ggplot2

<www.r-bloggers.com/basic-introduction-to-ggplot2/>

install.packages("ggplot2")
library(ggplot2)

Der diamonds Datensatz

head(diamonds)
carat cut color clarity depth table price x y z
0.23 Ideal E SI2 61.5 55 326 3.95 3.98 2.43
0.21 Premium E SI1 59.8 61 326 3.89 3.84 2.31
0.23 Good E VS1 56.9 65 327 4.05 4.07 2.31
0.29 Premium I VS2 62.4 58 334 4.20 4.23 2.63
0.31 Good J SI2 63.3 58 335 4.34 4.35 2.75
0.24 Very Good J VVS2 62.8 57 336 3.94 3.96 2.48

Wie nutzt man qplot

# histogram
qplot(depth, data=diamonds)

Ein Balkendiagramm

qplot(cut, depth, data=diamonds)

Ein weiteres Balkendiagramm

qplot(factor(cyl), data=mtcars,geom="bar")

Boxplot

qplot(data=diamonds,x=cut,y=depth,geom="boxplot")

Scatterplot

# scatterplot
qplot(carat, depth, data=diamonds)

Farbe hinzu:

qplot(carat, depth, data=diamonds,color=cut)

Trendlinie hinzufügen

myGG<-qplot(data=diamonds,x=carat,y=depth,color=carat) 
myGG + stat_smooth(method="lm")

Graphik drehen

qplot(factor(cyl), data=mtcars, geom="bar") + 
coord_flip()

Wie nutzt man ggplot

ggplot(diamonds, aes(clarity, fill=cut)) + geom_bar()

Farben selber wählen

Es wird das Paket RColorBrewer verwendet um die Farbpalette zu ändern

install.packages("RColorBrewer")
library(RColorBrewer)
myColors <- brewer.pal(5,"Accent")
names(myColors) <- levels(diamonds$cut)
colScale <- scale_colour_manual(name = "cut",
                                values = myColors)

http://stackoverflow.com/questions/6919025/

Eine Graphik mit den gewählten Farben

p <- ggplot(diamonds,aes(carat, depth,colour = cut)) + 
  geom_point()
p + colScale

Speichern mit ggsave

ggsave("Graphik.jpg")

Einfache Karten mit R erstellen

Gliederung

Arten von räumlichen Daten:

Das R-paket ggmap wird im folgenden genutzt um verschiedene Kartentypen darzustellen.

Mit qmap kann man eine schnelle Karte erzeugen.

Straßenkarten

Installieren des Paketes

devtools::install_github("dkahle/ggmap")
devtools::install_github("hadley/ggplot2")
install.packages("ggmap")

Paket ggmap - Hallo Welt

library(ggmap)

Und schon kann die erste Karte erstellt werden:

qmap("Mannheim")

Karte für eine Sehenswürdigkeit

BBT <- qmap("Berlin Brandenburger Tor")
BBT

Karte für einen ganzen Staat

qmap("Germany")

Ein anderes zoom level

qmap("Germany", zoom = 6)

Hilfe bekommen wir mit dem Fragezeichen

?qmap

Verschiedene Abschnitte in der Hilfe:

Die Beispiele in der Hilfe

Ausschnitt aus der Hilfe Seite zum Befehl qmap:

qmap Example

qmap Example

Das Beispiel kann man direkt in die Konsole kopieren:

# qmap("baylor university")
qmap("baylor university", zoom = 14)
# und so weiter

Ein anderes zoom level

qmap("Mannheim", zoom = 12)

Näher rankommen

qmap('Mannheim', zoom = 13)

Ganz nah dran

qmap('Mannheim', zoom = 20)

ggmap - maptype satellite

qmap('Mannheim', zoom = 14, maptype="satellite")

ggmap - maptype satellite zoom 20

qmap('Mannheim', zoom = 20, maptype="hybrid")

ggmap - maptype hybrid

qmap("Mannheim", zoom = 14, maptype="hybrid")

Terrain/physical maps

ggmap - terrain map

qmap('Schriesheim', zoom = 14,maptype="terrain")

Abstrahierte Karten (http://www.designfaves.com)

New York

New York

ggmap - maptype watercolor

qmap('Mannheim', zoom = 14,maptype="watercolor",source="stamen")

ggmap - source stamen

qmap('Mannheim', zoom = 14,
 maptype="toner",source="stamen")

ggmap - maptype toner-lite

qmap('Mannheim', zoom = 14,
 maptype="toner-lite",source="stamen")

ggmap - maptype toner-hybrid

qmap('Mannheim', zoom = 14,
 maptype="toner-hybrid",source="stamen")

ggmap - maptype terrain-lines

qmap('Mannheim', zoom = 14,
 maptype="terrain-lines",source="stamen")

Graphiken speichern

RstudioExport

RstudioExport

ggmap - ein Objekt erzeugen

MA_map <- qmap('Mannheim', 
               zoom = 14,
               maptype="toner",
               source="stamen")

Geokodierung

Geocoding (…) uses a description of a location, most typically a postal address or place name, to find geographic coordinates from spatial reference data …

Wikipedia - Geocoding

library(ggmap)
geocode("Mannheim",source="google")
lon lat
8.463243 49.48604

Latitude und Longitude

LatLon

LatLon

http://modernsurvivalblog.com

Koordinaten verschiedener Orte in Deutschland

cities lon lat
Hamburg 9.993682 53.55108
Koeln 6.960279 50.93753
Dresden 13.737262 51.05041
Muenchen 11.581981 48.13513

Reverse Geokodierung

Reverse geocoding is the process of back (reverse) coding of a point location (latitude, longitude) to a readable address or place name. This permits the identification of nearby street addresses, places, and/or areal subdivisions such as neighbourhoods, county, state, or country.

Quelle: Wikipedia

revgeocode(c(48,8))
## [1] "Unnamed Road, Somalia"

Die Distanz zwischen zwei Punkten

mapdist("Q1, 4 Mannheim","B2, 1 Mannheim")
##             from             to   m    km     miles seconds  minutes
## 1 Q1, 4 Mannheim B2, 1 Mannheim 749 0.749 0.4654286     212 3.533333
##        hours
## 1 0.05888889
mapdist("Q1, 4 Mannheim","B2, 1 Mannheim",mode="walking")
##             from             to   m    km     miles seconds minutes  hours
## 1 Q1, 4 Mannheim B2, 1 Mannheim 546 0.546 0.3392844     423    7.05 0.1175

Eine andere Distanz bekommen

mapdist("Q1, 4 Mannheim","B2, 1 Mannheim",mode="bicycling")
##             from             to   m    km    miles seconds  minutes
## 1 Q1, 4 Mannheim B2, 1 Mannheim 555 0.555 0.344877     215 3.583333
##        hours
## 1 0.05972222

Geokodierung - verschiedene Punkte von Interesse

POI1 <- geocode("B2, 1 Mannheim",source="google")
POI2 <- geocode("Hbf Mannheim",source="google")
POI3 <- geocode("Mannheim, Friedrichsplatz",source="google")
ListPOI <-rbind(POI1,POI2,POI3)
POI1;POI2;POI3
##        lon      lat
## 1 8.462844 49.48569
##        lon      lat
## 1 8.469879 49.47972
##        lon      lat
## 1 8.475208 49.48326

Punkte in der Karte

MA_map +
geom_point(aes(x = lon, y = lat),
data = ListPOI)

Punkte in der Karte

MA_map +
geom_point(aes(x = lon, y = lat),col="red",
data = ListPOI)

ggmap - verschiedene Farben

ListPOI$color <- c("A","B","C")
MA_map +
geom_point(aes(x = lon, y = lat,col=color),
data = ListPOI)

ggmap - größere Punkte

ListPOI$size <- c(10,20,30)
MA_map +
geom_point(aes(x = lon, y = lat,col=color,size=size),
data = ListPOI)

Eine Route von Google maps bekommen

from <- "Mannheim Hbf"
to <- "Mannheim B2 , 1"
route_df <- route(from, to, structure = "route")

Mehr Information

http://rpackages.ianhowson.com/cran/ggmap/man/route.html

Eine Karte mit dieser Information zeichnen

qmap("Mannheim Hbf", zoom = 14) +
  geom_path(
    aes(x = lon, y = lat),  colour = "red", size = 1.5,
    data = route_df, lineend = "round"
  )

Wie fügt man Punkte hinzu

http://i.stack.imgur.com

pic

pic

Cheatsheet

https://www.rstudio.com/

Cheatsheet

Cheatsheet

Resourcen und Literatur

Take Home Message

Was klar sein sollte:

Regressionsdiagnostik mit R-Paket visreg

Regressionsdiagnostik mit Basis-R

Ein einfaches Modell

N <- 5
x1 <- rnorm(N)
y <- runif(N)

Modellvorhersage machen

mod1 <- lm(y~x1)
pre <- predict(mod1)

Regressionsdiagnostik mit Basis-R

plot(x1,y)
abline(mod1)
segments(x1, y, x1, pre, col="red")

Das visreg-Paket

Ein Modell wird auf dem airquality Datensatz geschätzt

install.packages("visreg")
library(visreg)
fit <- lm(Ozone ~ Solar.R + Wind + Temp, data = airquality)

Visualisierung

visreg(fit)

Und dann mit visreg visualisiert.

visreg(fit, "Wind", type = "contrast")

Visualisierung mit dem Paket visreg

visreg(fit, "Wind", type = "contrast")

Das visreg-Paket

visreg(fit, "Wind", type = "conditional")

Regression mit Faktoren

Mit visreg können die Effekte bei Faktoren visualisiert werden.

airquality$Heat <- cut(airquality$Temp, 3, 
    labels=c("Cool", "Mild", "Hot"))
fit.heat <- lm(Ozone ~ Solar.R + Wind + Heat, 
    data = airquality)

Effekte von Faktoren

par(mfrow=c(1,2))
visreg(fit.heat, "Heat", type = "contrast")
visreg(fit.heat, "Heat", type = "conditional")

Das Paket visreg - Interaktionen

airquality$Heat <- cut(airquality$Temp, 3, 
labels=c("Cool", "Mild", "Hot"))
fit <- lm(Ozone ~ Solar.R + Wind * Heat, data = airquality)

Steuern der Graphikausgabe mittels layout

visreg(fit, "Wind", by = "Heat",layout=c(3,1))

Das Paket visreg - Interaktionen overlay

fit <- lm(Ozone ~ Solar.R + Wind * Heat, data = airquality)
visreg(fit, "Wind", by="Heat", overlay=TRUE, partial=FALSE)

Das Paket visreg - visreg2d

fit2 <- lm(Ozone ~ Solar.R + Wind * Temp, data = airquality)
visreg2d(fit2, "Wind", "Temp", plot.type = "image")

Das Paket visreg - surface

visreg2d(fit2, "Wind", "Temp", plot.type = "persp")

Weitere Themen im Ausblick

Mehr User Interface - Der R-commander

library("Rcmdr")

Interaktive Grafiken mit R

Das Paket iplots

install.packages("iplots",dep=TRUE)
library(iplots)
cyl.f <- factor(mtcars$cyl)
gear.f <- factor(mtcars$factor)
attach(mtcars)
ihist(mpg) # histogram
ibar(carb) # barchart
iplot(mpg, wt) # scatter plot
ibox(mtcars[c("qsec","disp","hp")]) # boxplots
ipcp(mtcars[c("mpg","wt","hp")]) # parallel coordinates
imosaic(cyl.f,gear.f) # mosaic plot

R-Paket rggobi

library(rggobi)
g <- ggobi(mydata)

Interaktion mit plots

attach(mydata)
plot(x, y) # scatterplot
identify(x, y, labels=row.names(mydata)) # identify points
coords <- locator(type="l") # add lines
coords # display list

Tabellen für Publikationen

library(stargazer)
stargazer(attitude)

Tabellen mit dem R-Paket knitr

library(knitr)
kable(head(iris), format = "latex")
pic

pic

Nichtlineare Regression

Folien zum Workshop:

https://github.com/Japhilko/npRegression/tree/master/slides

library(splines)

Beispiel einer interaktiven Karte

interaktiven Karte und Rcode um eine interaktive Karte mit leaflet zu erzeugen.